.ipynb .pdf

repository open issue

Binder

Contents

Manoeuvring simulations

Many simulation model for ship manoeuvring have been developed in the field of ship hydrodynamics such as: the Abkowitz model [Abk64] or the Norrbin model [Nor60]. This chapter will develop a general simulation model for ship manoeuvring, that can be further specified to become either the Abkowitz or Norbin model. Expressing the models on a general form is important in this research where many different models will be tested and compared.

The models can be expressed in a very general way as shown in eq_6DOF ( [Fos21]).

()\[\displaystyle C \nu + D \nu + M \dot{nu} + g + g_{0} = \tau + \tau_{wave} + \tau_{wind}\]

Where \(\eta\) describes the position and \(\nu\) is the velocities:

()\[\begin{split}\displaystyle \eta = \left[\begin{matrix}x_{0}\\y_{0}\\z_{0}\\\phi\\\theta\\\psi\end{matrix}\right]\end{split}\]

()\[\begin{split}\displaystyle \nu = \left[\begin{matrix}u\\v\\w\\p\\q\\r\end{matrix}\right]\end{split}\]

The model matrixes:

• \(M\) : inertia

• \(C\) : corriolis/centrepetal

• \(D\) : damping

• \(g(\eta)\) is a vector of generalized gravitational an buoyance forces.

• \(g_0\) is static restoring forces due to ballast systems.

The inertia as well as the corriolis matrix can be split into actual mass/intertia and added mass/inertia:

()\[\displaystyle M = M_{A} + M_{RB}\]

()\[\displaystyle C = C_{A} + C_{RB}\]

So that the model equation can be written as:

()\[\displaystyle C_{A} \nu_{r} + C_{RB} \nu + D \nu_{r} + M_{A} \dot{nu_r} + M_{RB} \dot{nu} + g + g_{0} = \tau + \tau_{wave} + \tau_{wind}\]

The velocity of the current can be described as (assumed to be irrotational):

()\[\begin{split}\displaystyle \nu_{r} = \left[\begin{matrix}u - u_{c}\\v - v_{c}\\w - w_{c}\\p\\q\\r\end{matrix}\right]\end{split}\]

So that \(\nu_r\) are the relative velocities (including current):

()\[\begin{split}\displaystyle \nu_{r} = \left[\begin{matrix}u - u_{c}\\v - v_{c}\\w - w_{c}\\p\\q\\r\end{matrix}\right]\end{split}\]

If the current is also assumed to be constant, acceleration of the relative velocities is the same:

()\[\displaystyle \dot{nu_r} = \dot{nu}\]

So that the entire equation can be expressed in terms of relative velocities:

()\[\displaystyle C_{A} \nu_{r} + C_{RB} \nu_{r} + D \nu_{r} + M_{A} \dot{nu_r} + M_{RB} \dot{nu_r} + g + g_{0} = \tau + \tau_{wave} + \tau_{wind}\]

Ship parameters

The Ocean Bird Wind Powered Car Carrier will be used as the reference ship of this book:

A model scale version of this ship has the following main dimensions:

{'T': 0.2063106796116504,
 'L': 5.014563106796117,
 'CB': 0.45034232324249973,
 'B': 0.9466019417475728,
 'rho': 1000,
 'x_G': 0,
 'm': 441.0267843660858,
 'I_z': 693.124396594905,
 'volume': 0.4410267843660858}

3 DOF

A simulation model with only three degrees of fredom (DOF) : surge, sway and yaw is often sufficient to describe the ship’s manoeuvring performance. The equation of motion can be expressed as (Triantafyllou & Hover, 2003):

\[\displaystyle m \left(\dot{u} - r^{2} x_{G} - r v\right) = X_{force}\]

\[\displaystyle m \left(\dot{r} x_{G} + \dot{v} + r u\right) = Y_{force}\]

\[\displaystyle I_{z} \dot{r} + m x_{G} \left(\dot{v} + r u\right) = N_{force}\]

Where \(X_{force}\), \(Y_{force}\) and \(N_{force}\) are the hydrodynamic forces from the ship. So these equations holds for most of the different models for manoeuvring, except linear models where the nonlinear \(r\cdot u\) term have been linearized. The difference in the methods is rather in how these hydrodynamic forces are calculated.

The hydrodynamic forces can be split into added masses (that depend on the accelerations: \(\dot{u}\), \(\dot{v}\), \(\dot{r}\)) and dampings (that depend on the velocities: \(u\), \(v\), \(r\)):

\[\displaystyle X_{force} = X_{\dot{u}} \dot{u} + \operatorname{X_{D}}{\left(u,v,r,\delta \right)}\]

\[\displaystyle Y_{force} = Y_{\dot{r}} \dot{r} + Y_{\dot{v}} \dot{v} + \operatorname{Y_{D}}{\left(u,v,r,\delta \right)}\]

\[\displaystyle N_{force} = N_{\dot{r}} \dot{r} + N_{\dot{v}} \dot{v} + \operatorname{N_{D}}{\left(u,v,r,\delta \right)}\]

\(X_{\dot{u}}\), \(X_{\dot{v}}\), \(X_{\dot{r}}\) are the added masses: when the ship is accelerating in water, not only the mass of the ship is accelerating, but also some mass in the water needs to be accelerated. This is referred to as added masses. Sometimes added masses such as \(X_{\dot{v}}\) is neglected, as in the case above.

The specialization of the various simulation models can now be further categorized to the functions to calculate the damping forces: \(X_{D}\), \(Y_{D}\), \(N_{D}\).

The general equation for all of the simulation models in this research can be written as:

\[\displaystyle m \left(\dot{u} - r^{2} x_{G} - r v\right) = X_{\dot{u}} \dot{u} + \operatorname{X_{D}}{\left(u,v,r,\delta \right)}\]

\[\displaystyle m \left(\dot{r} x_{G} + \dot{v} + r u\right) = Y_{\dot{r}} \dot{r} + Y_{\dot{v}} \dot{v} + \operatorname{Y_{D}}{\left(u,v,r,\delta \right)}\]

\[\displaystyle I_{z} \dot{r} + m x_{G} \left(\dot{v} + r u\right) = N_{\dot{r}} \dot{r} + N_{\dot{v}} \dot{v} + \operatorname{N_{D}}{\left(u,v,r,\delta \right)}\]

The classic simulation methods express the damping forces as a truncated taylor series. Here is a rather simple model, that has been created by the author to be used in the followin examples.

\[\displaystyle \operatorname{X_{D}}{\left(u,v,r,\delta \right)} = X_{delta} \delta + X_{rr} r^{2} + X_{r} r + X_{u} u + X_{vr} r v + X_{v} v\]

\[\displaystyle \operatorname{Y_{D}}{\left(u,v,r,\delta \right)} = Y_{delta} \delta + Y_{r} r + Y_{ur} r u + Y_{u} u + Y_{v} v\]

\[\displaystyle \operatorname{N_{D}}{\left(u,v,r,\delta \right)} = N_{delta} \delta + N_{r} r + N_{ur} r u + N_{u} u + N_{v} v\]

Similitude

Prime system

The prime system uses ship length \(L\), water density \(\rho\) and total velocity \(U^2=u^2+v^2\) to express physical quantities in nondimensional form. The quantities have a prime symbol (‘) attached to them when the prime system is used. The quantities are made nondimensional with the prime system by dividing with a denominator according to the table below:

	denominator
length	L
volume	L**3
mass	0.5*L**3*rho
density	0.5*rho
inertia_moment	0.5*L**5*rho
time	L/U
area	L**2
angle	1
-	1
linear_velocity	U
angular_velocity	U/L
linear_acceleration	U**2/L
angular_acceleration	U**2/L**2
force	0.5*L**2*U**2*rho
moment	0.5*L**3*U**2*rho

The manoeuvring models are often expressed in the prime system so that a nondimensional force \(F'\) for instance can be expressed as a function of some coefficients \(C\) and nondimensional velocity \(v\):

\[\displaystyle F' = C v'\]

This may look as a model where the force depend on the transverse velocity \(v\) only, but this is actually not true. If this equation is converted converted back to SI units it is instead written:

\[\displaystyle F = 0.5 C L^{2} U \rho v\]

So it is now clear that the example force model above, also depends on the total velocity \(U\).

Brix parameters

The hydrodynamic derivatives can be estimated with semi empirical formulas according to (Brix, 1993).

\[\displaystyle X_{\dot{u}} = 0.00017880704448812\]

\[\displaystyle Y_{\dot{v}} = -0.00610938740826336\]

\[\displaystyle Y_{\dot{r}} = -0.000303137774581419\]

\[\displaystyle N_{\dot{v}} = -0.000128254401723757\]

\[\displaystyle N_{\dot{r}} = -0.000298674827731788\]

\[\displaystyle X_{u} = -0.001\]

\[\displaystyle X_{v} = 0.0\]

\[\displaystyle X_{r} = 0.0\]

\[\displaystyle X_{delta} = -0.001\]

\[\displaystyle X_{vr} = -0.006\]

\[\displaystyle X_{rr} = 0.007\]

\[\displaystyle Y_{u} = 0.0\]

\[\displaystyle Y_{v} = -0.00971290834761524\]

\[\displaystyle Y_{r} = 0.00240236113642936\]

\[\displaystyle Y_{delta} = 0.001\]

\[\displaystyle Y_{ur} = 0.001\]

\[\displaystyle N_{u} = 0.0\]

\[\displaystyle N_{v} = -0.00318395170922333\]

\[\displaystyle N_{r} = -0.00171885123535641\]

\[\displaystyle N_{delta} = -0.0005\]

\[\displaystyle N_{ur} = 0.0\]

The figure above shows the quasi static forces for the present ship with the hydrodynamic derivatives during a variation of rudder angle \(\delta\), yaw rate \(r\) and tranverse velocity \(v\) and total velocity \(U\)=2 m/s

Simulation

Decoupling

There is a coupling between the sway and yaw equation thruogh the added masses. These equations need to be decoupled to be usable in a simulation model. This is done by first expressing the equation in matrix form:

\[\begin{split}\displaystyle \left[\begin{matrix}- X_{\dot{u}} + m & 0 & 0\\0 & - Y_{\dot{v}} + m & - Y_{\dot{r}} + m x_{G}\\0 & - N_{\dot{v}} + m x_{G} & I_{z} - N_{\dot{r}}\end{matrix}\right] \left[\begin{matrix}\dot{u}\\\dot{v}\\\dot{r}\end{matrix}\right] = \left[\begin{matrix}- m \left(- r^{2} x_{G} - r v\right) + \operatorname{X_{D}}{\left(u,v,r,\delta \right)}\\- m r u + \operatorname{Y_{D}}{\left(u,v,r,\delta \right)}\\- m r u x_{G} + \operatorname{N_{D}}{\left(u,v,r,\delta \right)}\end{matrix}\right]\end{split}\]

The decoupled equations are obtained by inverting the intertia matrix:

This equation can be written in a cleaner way if the \(S\) term is introduced:

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{u}\\\dot{v}\\\dot{r}\end{matrix}\right] = \left[\begin{matrix}\frac{1}{- X_{\dot{u}} + m} & 0 & 0\\0 & - \frac{- I_{z} + N_{\dot{r}}}{S} & - \frac{- Y_{\dot{r}} + m x_{G}}{S}\\0 & - \frac{- N_{\dot{v}} + m x_{G}}{S} & - \frac{Y_{\dot{v}} - m}{S}\end{matrix}\right] \left[\begin{matrix}- m \left(- r^{2} x_{G} - r v\right) + \operatorname{X_{D}}{\left(u,v,r,\delta \right)}\\- m r u + \operatorname{Y_{D}}{\left(u,v,r,\delta \right)}\\- m r u x_{G} + \operatorname{N_{D}}{\left(u,v,r,\delta \right)}\end{matrix}\right]\end{split}\]

Where \(S\) is a helper variable defined as:

\[\displaystyle S = - I_{z} Y_{\dot{v}} + I_{z} m + N_{\dot{r}} Y_{\dot{v}} - N_{\dot{r}} m - N_{\dot{v}} Y_{\dot{r}} + N_{\dot{v}} m x_{G} + Y_{\dot{r}} m x_{G} - m^{2} x_{G}^{2}\]

prime and SI

This concise expression to calculate the accelerations is expressed in the prime system. Conducting a simulation by integrating this expression will give a simulation that is conducted in the prime system, with prime system time etc., which is not very convenient. If a simulation should instead be conducted with a standard unit system (SI), the states (velocities as well as the resulting accelerations) and the ship parameters (length, masses and inertia) need to be converted.

The input velocities in SI should be converted to prime and the output accelerations in prime needs to be converted back to SI:

1. Input: velocities : SI -> prime

2. Output: accelerations : prime -> SI

This can be done for each time step in the simulations, by implementing this logic into the code. Another alternative is to rewrite the expression into SI units, which will be done below. The damping functions: \(X_D(u,v,r,\delta)\), \(Y_D(u,v,r,\delta)\), \(N_D(u,v,r,\delta)\) is also substituted with their polynomials to derive the full expression.

u1d_function = sp.Function(r'\dot{u}')(u,v,r,delta)
v1d_function = sp.Function(r'\dot{v}')(u,v,r,delta)
r_function = sp.Function(r'\dot{r}')(u,v,r,delta)

subs_prime = [
    
    (m,m/prime_system.df_prime.mass.denominator),
    (I_z,I_z/prime_system.df_prime.inertia_moment.denominator),
    (x_G,x_G/prime_system.df_prime.length.denominator),
    
    (u, u/sp.sqrt(u**2+v**2)),
    (v, v/sp.sqrt(u**2+v**2)),
    (r, r/(sp.sqrt(u**2+v**2)/L)),
   
]

subs = [
(X_D, vmm.X_qs_eq.rhs),
(Y_D, vmm.Y_qs_eq.rhs),
(N_D, vmm.N_qs_eq.rhs),
]

subs = subs + subs_prime

A_SI = A.subs(subs)
b_SI = b.subs(subs)

x_dot = sympy.matrices.dense.matrix_multiply_elementwise(A_SI.inv()*b_SI,
                                                 sp.Matrix([(u**2+v**2)/L,(u**2+v**2)/L,(u**2+v**2)/(L**2)]))

sp.Eq(sp.UnevaluatedExpr(sp.Matrix([
    u1d_function,
    v1d_function,
    r_function,    
])),sp.UnevaluatedExpr(x_dot))

\[\begin{split}\displaystyle \left[\begin{matrix}\dot{u}{\left(u,v,r,\delta \right)}\\\dot{v}{\left(u,v,r,\delta \right)}\\\dot{r}{\left(u,v,r,\delta \right)}\end{matrix}\right] = \left[\begin{matrix}\frac{1.0 L^{2} \rho \left(u^{2} + v^{2}\right) \left(\frac{L^{2} X_{rr} r^{2}}{u^{2} + v^{2}} + \frac{L X_{r} r}{\sqrt{u^{2} + v^{2}}} + \frac{L X_{vr} r v}{u^{2} + v^{2}} + X_{delta} \delta + \frac{X_{u} u}{\sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + \frac{X_{v} v}{\sqrt{u^{2} + v^{2}}} - \frac{2.0 m \left(- \frac{L r^{2} x_{G}}{u^{2} + v^{2}} - \frac{L r v}{u^{2} + v^{2}}\right)}{L^{3} \rho}\right)}{- 1.0 L^{3} X_{\dot{u}} \rho + 2.0 m}\\\frac{\left(u^{2} + v^{2}\right) \left(\frac{\left(- 4.0 I_{z} L^{3} \rho + 2.0 L^{8} N_{\dot{r}} \rho^{2}\right) \left(\frac{L Y_{r} r}{\sqrt{u^{2} + v^{2}}} + \frac{L Y_{ur} r u}{\sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + Y_{delta} \delta + \frac{Y_{u} u}{\sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + \frac{Y_{v} v}{\sqrt{u^{2} + v^{2}}} - \frac{2.0 m r u}{L^{2} \rho \sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}}\right)}{4.0 I_{z} L^{3} Y_{\dot{v}} \rho - 8.0 I_{z} m - 2.0 L^{8} N_{\dot{r}} Y_{\dot{v}} \rho^{2} + 2.0 L^{8} N_{\dot{v}} Y_{\dot{r}} \rho^{2} + 4.0 L^{5} N_{\dot{r}} m \rho - 4.0 L^{4} N_{\dot{v}} m \rho x_{G} - 4.0 L^{4} Y_{\dot{r}} m \rho x_{G} + 8.0 m^{2} x_{G}^{2}} + \frac{\left(- 2.0 L^{8} Y_{\dot{r}} \rho^{2} + 4.0 L^{4} m \rho x_{G}\right) \left(\frac{L N_{r} r}{\sqrt{u^{2} + v^{2}}} + \frac{L N_{ur} r u}{\sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + N_{delta} \delta + \frac{N_{u} u}{\sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + \frac{N_{v} v}{\sqrt{u^{2} + v^{2}}} - \frac{2.0 m r u x_{G}}{L^{3} \rho \sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}}\right)}{4.0 I_{z} L^{3} Y_{\dot{v}} \rho - 8.0 I_{z} m - 2.0 L^{8} N_{\dot{r}} Y_{\dot{v}} \rho^{2} + 2.0 L^{8} N_{\dot{v}} Y_{\dot{r}} \rho^{2} + 4.0 L^{5} N_{\dot{r}} m \rho - 4.0 L^{4} N_{\dot{v}} m \rho x_{G} - 4.0 L^{4} Y_{\dot{r}} m \rho x_{G} + 8.0 m^{2} x_{G}^{2}}\right)}{L}\\\frac{\left(u^{2} + v^{2}\right) \left(\frac{\left(- 1.0 L^{8} N_{\dot{v}} \rho^{2} + 2.0 L^{4} m \rho x_{G}\right) \left(\frac{L Y_{r} r}{\sqrt{u^{2} + v^{2}}} + \frac{L Y_{ur} r u}{\sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + Y_{delta} \delta + \frac{Y_{u} u}{\sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + \frac{Y_{v} v}{\sqrt{u^{2} + v^{2}}} - \frac{2.0 m r u}{L^{2} \rho \sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}}\right)}{2.0 I_{z} L^{3} Y_{\dot{v}} \rho - 4.0 I_{z} m - 1.0 L^{8} N_{\dot{r}} Y_{\dot{v}} \rho^{2} + 1.0 L^{8} N_{\dot{v}} Y_{\dot{r}} \rho^{2} + 2.0 L^{5} N_{\dot{r}} m \rho - 2.0 L^{4} N_{\dot{v}} m \rho x_{G} - 2.0 L^{4} Y_{\dot{r}} m \rho x_{G} + 4.0 m^{2} x_{G}^{2}} + \frac{\left(1.0 L^{8} Y_{\dot{v}} \rho^{2} - 2.0 L^{5} m \rho\right) \left(\frac{L N_{r} r}{\sqrt{u^{2} + v^{2}}} + \frac{L N_{ur} r u}{\sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + N_{delta} \delta + \frac{N_{u} u}{\sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}} + \frac{N_{v} v}{\sqrt{u^{2} + v^{2}}} - \frac{2.0 m r u x_{G}}{L^{3} \rho \sqrt{u^{2} + v^{2}} \sqrt{u^{2} + \frac{v^{2}}{u^{2} + v^{2}}}}\right)}{2.0 I_{z} L^{3} Y_{\dot{v}} \rho - 4.0 I_{z} m - 1.0 L^{8} N_{\dot{r}} Y_{\dot{v}} \rho^{2} + 1.0 L^{8} N_{\dot{v}} Y_{\dot{r}} \rho^{2} + 2.0 L^{5} N_{\dot{r}} m \rho - 2.0 L^{4} N_{\dot{v}} m \rho x_{G} - 2.0 L^{4} Y_{\dot{r}} m \rho x_{G} + 4.0 m^{2} x_{G}^{2}}\right)}{L^{2}}\end{matrix}\right]\end{split}\]

…so the expanded expression in SI units is a lot more complicated.

Simulate data

\[\displaystyle \dot{\vec{x}} = f{\left(\vec{x},u_{input},w_{noise} \right)}\]

Where the state vector \(\vec{x}\):

\[\begin{split}\displaystyle \vec{x} = \left[\begin{matrix}x_{0}\\y_{0}\\\Psi\\u\\v\\r\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle f{\left(\vec{x},u_{input},w_{noise} \right)} = \left[\begin{matrix}u \cos{\left(\Psi \right)} - v \sin{\left(\Psi \right)}\\u \sin{\left(\Psi \right)} + v \cos{\left(\Psi \right)}\\r\\\dot{u}{\left(u,v,r,\delta \right)}\\\dot{v}{\left(u,v,r,\delta \right)}\\\dot{r}{\left(u,v,r,\delta \right)}\end{matrix}\right]\end{split}\]

The actual simulation can be carried out with a simple Euler forward integration:

def time_step(x_,u_):
    psi=x_[2]
    u=x_[3]
    v=x_[4]
    r=x_[5]
    delta = u_
    x_dot = run(lambda_f, **parameters, **ship_parameters, psi=psi, u=u, v=v, r=r, delta=delta).flatten()
    return x_dot

def simulate(x0,E, ws, t, us):
    
    simdata = np.zeros((6,len(t)))
    x_=x0
            
    for i,(u_,w_) in enumerate(zip(us,ws)):
       
        x_dot = time_step(x_,u_)
        
        x_=x_ + h_*x_dot
        
        simdata[:,i] = x_
                
    df = pd.DataFrame(simdata.T, columns=["x0","y0","psi","u","v","r"], index=t)
    df.index.name = 'time'
    df['delta'] = us
    
    return df

parameters=df_parameters['prime'].copy()
N_ = 1000
t_ = np.linspace(0,50,N_)
h_ = float(t_[1]-t_[0])
ws = np.zeros(N_)
us = delta_ = np.deg2rad(20)*np.ones(N_)
x0_ = np.array([0,0,0,2,0,0])

E=np.array([0,0,0,1,1,1])
df = simulate(x0=x0_, E=E, ws=ws, t=t_, us=us)
df['beta'] = -np.arctan2(df['v'],df['u'])
df['U'] = np.sqrt(df['u']**2 + df['v']**2)
df.to_csv('first_simulation.csv', index=True)

Simulate parameter contributions

Here is an interavtive graph showing how the various hydrodynamic derivatives contribute to the forces:

previous

Introduction example

next

Manoeuvring parameter identification